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Abstract 

We consider a direct optimization approach for ensemble density func- 
tional theory electronic structure calculations. The update operator for 
the electronic orbitals takes the structure of the Stiefel manifold into ac- 
count and we present an optimization scheme for the occupation numbers 
that ensures that the constraints remain satisfied. Wc also compare se- 
quential and simultaneous quasi-Newton and nonlinear conjugate gradient 
optimization procedures, and demonstrate that simultaneous optimization 
of the electronic orbitals and occupation numbers improve performance 
compared to the sequential approach. 

1 Introduction 

Advances in computer power and numerical methods during the past few decades 
has dramatically increased the scope of electronic structure problems that can be 
computationally studied. Kohn-Sham density functional theory (DFT) meth- 
ods can be used to reach precision comparable to experimental accuracy for 
insulators and semiconductors, while metallic systems remain more challeng- 
ing. Metallic systems lack a gap between occupied and unoccupied electronic 
states in the energy spectrum, which leads to slower convergence compared to 
insulators and semiconductors. Smearing of the Fermi surface is often used to 
enable convergence of metallic systems as well as insulators at positive temper- 
atures. Ensemble DFT permits direct computation of the occupation numbers 
of the orbitals based on the entropic term in the Helmholtz free energy. We 
consider an optimization problem where the target functional A corresponds to 
the Helmholtz free energy and the variables X and f to the electronic orbitals 
and occupation numbers respectively. 
The optimization problem is therefore 

minimize A(X, f), (1) 

* Department of Mathematics and Systems Analysis, Aalto University School of Science, 
Espoo, Finland, e-mail: kurt.baarman9aalto.fi 

t Department of Applied Physics, Aalto University School of Science, Espoo, Finland 



subject to 

X e = {X e M"^" I X'^X = I}. (2) 

Furthermore f G M" with fi ~ ^^'^ < /i < 1, where ne G N is the 

number of electrons in the system and Ue < n. We also assume that Vx^(X, f ) 
and Vf A(X, f ) are available, but expensive to compute. However, due to the 
form of A(X, f ) the price to compute A, Vx^, and Vf A simultaneously is 
comparable to computing one of them separately. Furthermore we assume that 
m ^ n, and that m is sufficiently large as to make storage of and operation 
with full m X m matrices prohibitively expensive. 

The orthogonality constraint on X means that A4 C M™^" defines the Stiefel 
manifold, which has the tangent space 

rxAl = {Y = XA + Z| A'^ = -AandZ'^X = 0}, (3) 

where Y, Z e " and A e M"^". We use the standard inner product 

(X,Y) =trace(X'^Y). (4) 

Given an arbitrary matrix W £ jj^x" -^ve can orthogonally project it onto Tx-A^ 
with 

Y = Px(W) = (I- iXX'^)W- iXW^X. (5) 

Minimization approaches to non-temperature dependent DFT do not in gen- 
eral permit fractional occupation of electronic orbitals [4| [l5|p4] - [26|[28] . In con- 
trast, explicit minimization with regards to occupation numbers permits frac- 
tional occupation based on the entropy functional of the Helmholtz free energy 
and can improve convergence, especially for metallic systems 6 p3p9 . It is also 



possible to transform Equation ([!]) into a nonlinear eigenvalue problem that can 
be solved through a self consistent field iteration 15 18 23 24 . The absence 
of well separated occupied and unoccupied orbitals make metallic systems chal- 
lenging to compute, and broadening of the Fermi surface is used to facilitate 
convergence [5j|27]. This broadening is often achieved by assigning the orbitals 
close to the Fermi level a fractional occupation number determined by the en- 



ergy of the electronic orbital 14 20 21 . Direct minimization on the other hand 
does not require the orbital energies to be computed at every step, and these 
broadening schemes are therefore not well suited for minimization methods. 



In 11 a framework for optimization methods on the Stiefel and Grassmann 
manifolds is presented, while [oj discusses a Newton-like iteration scheme on a 
more general manifold. Univariate optimization methods for the Stiefel manifold 
is presented in |7| , where identity plus rank one Householder transforms are given 
as one possible choice for moving on the manifold. The choice of coordinates 
can also be based on a QR factorization and polar decompositions [8yl0) or Lie 



groups 16 . An overview of geometric numerical integration techniques can be 
found in |17|. 

In Section [2] we first recall the nonlinear conjugate gradient and the quasi- 
Newton methods adapted for use on the Stiefel manifold. We then present 
an optimization procedure for the occupation numbers and end the section by 
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presenting a simultaneous orbital-occupation optimization strategies. Then, in 
Section [3] we numerically demonstrate the method on a model problem that 
includes nonlincarities similar to a DFT problem. The conclusions are finally 
presented in Section [4j 



2 Optimization with orthogonality constraints 
2.1 Update and transport 

We ensure that X^+i satisfies the orthogonality constraint by using a unitary 
update operator U which maps — > 7M. A search direction Y e Tx-^ given 
by an optimization procedure can be written 



Y = XA + QR, 



(6) 



where Q e R*"^", A,R e K"''", A^ = -A, Q^Q = I, and Q^X = 0. If the 
terms in Equation (|6| are not full rank the size of the matrices can be adjusted 
accordingly. 

If we follow Y to update X along a Stiefel geodesic we obtain the update 
operator for X ^llj 



U = [X Q] exp ( T 



A 
R 



R^ 




(7) 



with step length parameter t. The update operator generalized for an arbitrary 
matrix in span(X, Q) is 
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[X Q]- 



(8) 



where the orthogonality of X and Q has been exploited. 

In order to use information gained from previous evaluations of A and 
we must take M. into account. This requires us to transport vectors Y S Tx-A^ 
to 7ux-A4 with the transport operator 



T = I,„ + [X Q] exp T 



A R^ 
R 



I 



2ri 



[X Q]' 



(9) 



Here I„, e 



and lon e 



p2nx2n 



, and T does not modify matrices Z that 



satisfy [X Q] Z = 0. 

Remark 1: The closely related Grassmann manifold is identical to the 
Stiefel manifold with the addition of the homogeneity condition A(X) = A(XQ), 
where Q is orthogonal. The homogeneity condition is satisfied for orbitals with 
identical occupation numbers, but does not generally hold for ensemble DFT. A 
discussion of direct minimization with integer occupation numbers is presented 

ing. 
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2.2 Nonlinear conjugate gradients 

The conjugate gradient (CG) method can be viewed as an optimization method 
for a quadratic problem. Several generalizations of the CG method have been 
presented to solve optimization problems that are not quadratic [22| . Below, we 
review a nonlinear CG (NLCG) method adapted to account for the curvature 



of the manifold 11 



Given Xq which satisfies Xq Xq = I, the gradient projected onto Txo-^ is 

Fo = Pxo(Vx^(Xo,fo)), (10) 

and the initial search direction is the direction of steepest descent 

Yo = -Fo. (11) 

On the manifold the NLCG method then proceeds by minimizing A along 
the path defined by the search direction Y^. In practice we evaluate A once 
along the search direction and construct a quadratic approximation that we 
minimize. The step length, rj,, that minimizes A along the search direction is 
then used to update X^ such that 

Xfe+i = T(rfc)Xfe, (12) 

and the gradient and search directions are transported to Tx^^^Ai by T(rfc). 
The new projected gradient 

Ffc+i = Px,+,(VxA(X,+i,ffc+i)), (13) 

and search direction 

Yfc+i ^ -Fk+i + 7fcT(rfe)Yfc, (14) 
are then computed where 

(Ffc+i - T(rfe)Ffe,Ffe+i) 

= • ^''^ 

The step length is determined by the minimizer of a quadratic approximation 
of A along the search direction. The quadratic approximation is constructed by 
taking a trial step length Tg = inax(T,nin, Tfe-i), where Tmin is a predefined 
minimum trial step length and computing 

p(0)=A(X,f), 

p(re)=A(T(re)X,f), (16) 
p'(0) = (Y,Vx^(X,f)). 

Then solve Tk and limit it by 2rfc_i, and construct the update T(Tfc). This 
approximate line search requires one extra evaluation of A per step. 
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2.3 Quasi-Newton method 

The quasi-Newton (QN) method is similar to Newton's method, but replaces 
the inverse Hessian with an approximation. This is frequently possible even 
when the Hessian is not available, and can still be used to improve performance 
for a badly conditioned minimization problem. 

We base the QN method on Broyden's second or bad generalized update to 
construct the approximate inverse Hessian, G, of A at X^. While Broyden's 
second update does not construct a symmetric approximation, or ensure that 
the approximation is positive definite it is a robust update choice for electronic 
structure calculations (2[[l8]. Furthermore, X and Vx^ are matrices, 
which we take into account when constructing the generalized Broyden update. 
The secant condition is then 

GA* = AH, (17) 

where A$ and AH are the collected orbital gradient and position differences 
projected onto the tangent space and transported to Tk^AA. That is 

A*=[AFfe_i T(Tfc_i)AF,_2 ... T(T,_i)...T(ri+i)AF,] , (18) 

and 

AH=[AXfc_i T(Tfc_i)AXfc_2 ... T(rfe_i)...T(Tz+i)AX,] , (19) 

for history length k — I. Here the gradient differences projected onto Txt+i-M 
are 

AF, - F,+i - T(rOF„ (20) 

and Fi is like in (131, 

F, =Px,(Vx^(X„f,)). (21) 
The projected occupation weighted orbital differences are 

AX, - Px.+, (X,+i diag(f,+i) - X, diag(f,)) , (22) 

and the motivation for including the weight is that the unoccupied electronic 
orbitals do not contribute to the energy of the system. The no change condition 
is now 

Z = GZ VZ such that Z^A* = 0. (23) 

The secant and no change condition together correspond to the generalized Broy- 
den's second update where all single orbital secant conditions are simultaneously 
enforced for the entire history length. We can therefore use the generalized up- 



date formula 12 



G = ^1+ (AH-/xA*)(A*^A$)-iA*^, (24) 

where dropping the empty orbitals ensure that A*^A* is nonsingular in prac- 
tice. The search direction given by the QN method is then 

Y = -GF, (25) 
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and 

Xfe+i =U(rfc)Xfc, (26) 

where Y determines U as in Section ( |2.1| . The hne search is identical to the 
one described for the NLCG method in Section 12.21 with the addition of the 
constant underrelaxation /3x G ]0: 1] that we have included in the step length 

Tfc. 

In practice only the last few history steps contribute significantly to the rate 
of convergence. Consequently, we discard the oldest trial solutions and gradient 
information when a predetermined history length is reached. 



2.4 Optimization of occupation numbers 

Given a set of electronic orbitals X it is possible to further reduce A by optimiz- 
ing f . Forcing occupation towards a uniform distribution increases contributions 
to A from higher energy states, while simultaneously increasing the entropy 
which contributes to a reduction of A at nonzero temperatures. The relative 
strength of both of these effects determine the ground state of the system, and 
can lead to nonzero occupation of higher energy states at positive temperatures 
or due to nonlinear effects. 

Therefore, given X, we want to find f that minimizes A. To keep the number 
of particles constant we determine the search direction y which is the vector 
closest — Vf A(X, f ) that ensures that the conditions J2^i=i fi — "-e and < fi < 
1 remain satisfied. To this end we solve 



minimize ||y + VfA(X,f)||, (27) 

with the constraints X]"=i Vi — Hi < if /i = 1, and yi > if /i = 0. The first 
constraint on y ensures that the minimization step conserves electrons while the 
second and third condition prohibits unphysical occupation numbers. In prac- 
tice we use the quadprog routine available in MATLAB to solve this problem. 
Given the search direction y we minimize A by constructing a quadratic ap- 



proximation similar to ( 16 ) 



After we have solved y the occupation step length is determined like in 



Section 2.2 with the addition of the constant underrelaxation /3f S ]0, 1] included 
in (Tfc. In addition, we ensure that the occupation remains physical by limiting 
ak with ctm such that < fi + auVi < 1 for all i.lt is possible to take a longer 



step than ctm by recomputing y from Equation ( 27 1 with the updated boundary 
information when an entry in f reaches the boundary of physical occupation, 
or 1. However, convergence of occupation numbers is faster than orbital conver- 
gence, and the numbers of steps needed for convergence is therefore determined 
by the orbital convergence. Furthermore, if the occupation number of the least 
populated orbital has been less than 10^^^ on two consecutive iterations we 
drop the associated orbitals. 
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2.5 Simultaneous step size selection 

Typically an ensemble DFT problem is solved by sequentially optimizing the or- 
bitals with fixed occupation numbers and then fixing the orbitals and optimizing 
the occupation numbers. This process is then repeated until a satisfactory so- 
lution is obtained. 

The cost of evaluating A, Vx^, and Vf A is comparable to evaluating one 
of them separately, and simultaneous optimization of A with respect to X and 
f can for this reason potentially reduce computational effort. 

Given a pair of search directions (Y, y) for the orbitals and occupation 
numbers respectively and starting guesses for step lengths, Tk-i and ak-i we 
evaluate A and its gradients with the following trial step lengths 

Te = jQ niax(rmin, Tk-i), (28) 

and 

o-g = min(crM, 3^ max(CT„iin, crfc_i)). (29) 

Here Tmin and (Tmin are minimum trial step lengths. With this we construct a 
quadratic surface approximation 

p{t, a) = CiT^ + C2cr^ + c^t + 04(7 + C5, (30) 

that we use to simultaneously update both X and f by evaluation in one trial 
point. This surface is determined by the system of equations 

p(0,0)=v4(X,f), 
p,(0,0) = (Y,VxA(X,f)), 

P.(0,0) = (yo,VfA(X,f)), (31) 

PriTe, Ue) = (Y, Vxv4(T(re)X, f + G,y)), 
Pa{Te,<Je) = {ya^,WfA{T{T,)X,i + <J,y)). 

Solving this system and finding the minimums gives the optimal step lengths 
ffc and (Tfc for the quadratic approximation of the search directions. For the 
simultaneous NLCG method the step lengths are then Tk = and ak = 
while the QN method uses Tk = l3xTk and ak = Pf<Jk, where /3x,/3f £ ]0, 1] are 
constant underrelaxation parameters. We then simultaneously update X and f 
with Xfc+i = T(Tfc)Xfc and ffc+i = ffc + min(crM, o-fc)y respectively. 

We then simultaneously update X and f with Xfc+i = T(rfc)Xfc and ffc+i = 
ffc + min((TM, o-fc)y respectively. 



Remark 2: The surface (30 1 is determined by computing Vx^ and VfA at 



the trial step. The system of equations (31) could alternatively be determined 
by computing both A and Vx^ or A and VfA at (T(Te)X,f + <Jey). 

Remark 3: Inclusion of the tct cross term would require an extra trial 



evaluation point for system (31) to be linearly independent 
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3 Numerical experiments 



We use a two dimensional model problem to compare the sequential and simul- 
taneous NLCG and QN methods. This model problem is inspired by ensemble 
DFT, and corresponds to a three dimensional system constrained to two dimen- 
sions without spin effects and exchange-correlation terms while taking entropy 
into account. The model problem adapted from Reference [l9] is 

^(X,f) = -itr(X^LXdiag(f)) -t-vj,tn+ iv;r.tn-r5(f). (32) 

Here L G ^rnxm ^Yie discretized Laplace operator, Voxt G M™ the external 
potential, n e M™ the electron density, Vjnt — Vn the Hartree potential cor- 
responding to the electron density n, T to temperature, and S is the entropy. 
The electron density is 

n = (X o X) f , (33) 
where o is the entrywise, or Hadamard, product. The entropy term is 

n 

S{i) = - 5] /. ln(/, + <5(1 - /,)) + (1 - /,) ln(l - /, + Sf,), (34) 
1=1 

where (5 > is a small regularization parameter that ensures that the derivative 
of S remains finite. 

To calculate the potentials we use 



N 



(ve.t). = -V^^^V^' (35) 



where the sum is over the nuclei with charge Zj and position Rj. The position 
corresponding to the discretization point i is r^, and the parameter a is used to 
regularize the potential. V G M'"'^"' is similarly given by 

Vi," = T, K, • (36) 

||r, -r,||+a ^ ' 

We solve the problem in the unit square with zero boundary conditions 
corresponding to an infinite potential well. We use a uniform finite difference 
discretization with m inner points to obtain a system where X G E™^". Here 
n corresponds to the number of electronic orbitals in the calculation. As initial 
guess we use the solution of the quadratic problem using the first two terms 



of ( 32 1 . The occupation numbers are initialized to 

+ (37) 

where A — min(ne/n, 1 — rie/n) and Ue < n is the number of electrons. This 
choice ensures that the initial occupation of all orbitals is nonzero and empha- 
sizes lower energy orbitals. 
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We demonstrate the methods for three external potentials. For all models we 
use potential regularization a = 5 x 10~^ and entropy regularization 6 = 10"'^. 

The first model is a single nucleus with charge Z — 2 centered at the center of 
the unit square with two electrons. For this system the second and third orbitals 
are degenerate. We calculate the model with 10 electronic orbitals and a first 
order finite difference discretization with 25 interior points in one dimension 
resulting in rn = 625 spatial degrees of freedom. We will refer to this system as 

Z2. 

The second model, which we name Z3-Z2, consists of two nuclei, with a 
nuclei of charge Z = 3 placed at (|, |) and another with charge Z — 2 placed 
at |, I and 5 electrons. This system has four well separated electronic orbitals, 
while the fifth and sixth are relatively close. The computation is initialized with 
10 orbitals and 29 interior grid points in one dimension for m = 841. 

The last model, Z4-Z3, consists of two nuclei, Z = A placed at (|, |), and 
Z — 3 a.t the grid point closest to (|, ^) and 7 electrons. The off diagonal 
placement is chosen to break the symmetry of the system. This model initially 
has 14 orbitals and 29 interior grid points in one dimension (m = 841). 

For the sequential QN orbital minimizer uses the parameters /3x = 0.4, 
^ — 5 X 10~^, and history length 6. The sequential QN and NLCG methods 
minimum trial step length tq — 10~^ and we perform 6 orbital optimization 
steps before engaging the occupation number minimizer. Both sequential op- 
timization routines use an identical SD routine with /3f = 0.5, /i — lO"'* and 
(To = 10"'* for occupation number optimization with two optimization steps. We 
have tried several different combinations of orbital and occupation optimization 
steps and observed that this combination offers a good compromise. For the SD 
method /x only serves to scale the approximate line search. 

We measure convergence by the energy difference to a reference energy com- 
puted by running the simultaneous methods for 3000 steps and the sequential 
methods for 3000 optimization rounds. We then use the lowest energy obtained 
as the reference energy. 

The change in occupation numbers with rising temperature is graphically 
presented in Figures [T][2j and [3] The same data is repeated in Tables [TJ [2] and[3j 
At T = the lowest electronic orbitals are fully occupied for Z4-Z3, while one 
electron is split between two degenerate orbitals for Z2. Even though there is a 
small gap (1.69 x 10"^) between the fifth and sixth electron orbitals for Z3-Z2 
the fifth electron is split (0.55 vs 0.45) between these orbitals. We successfully 
replicated this split with a 3000 round sequential SD orbital occupation number 
optimization. Furthermore, restarting the SD iteration with a five orbital initial 
guess based on the split orbital reference solution results in convergence to a 
higher energy state. 

Figures |4j [5) and |6] illustrate energy convergence for the different methods. 
The simultaneous methods generally perform better than the sequential meth- 
ods, and the simultaneous NLCG method is more robust than the simultaneous 
QN approach. In the energy convergence for the sequential optimization rou- 
tines the switch between orbital and occupation optimization is readily seen in 
the steplike energy convergence. Furthermore, the performance of the sequential 
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Figure 1: Orbital energy levels with occupation for Z2 at varying temperatures. 
Fractional occupation numbers are indicated and the same data is also presented 
in Tabled] 
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Figure 2: Orbital energy levels with occupation for Z^-Z^ at varying temper- 
atures with orbital energy shifted by +5. Fractional occupation numbers are 
indicated and the same data is also presented in Table [2] 
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Figure 3: Orbital energy levels with occupation for at varying temper- 

atures with orbital energy shifted by +10. Fractional occupation numbers are 
indicated and the same data is also presented in Table |3] 



QN and NLCG methods is nearly identical for all models. This might be due to 
the limited number of step available for orbital optimization before occupation 
optimization is enabled. 

The simultaneous NLCG method outperforms the QN method for the Z4-Z3 
system shown in Figure [6j Increasing history generally improves the convergence 
rate of the QN method, but this did not significantly change the rate of con- 
vergence for this model. Frequent restarts limit history length and provide at 
least a partial explanation for this effect. Figure [7] presents Z4-Z3 restarts for 
the simultaneous QN method. Restarts are frequent for this model at all tem- 
peratures compared to Z2 and Z3-Z2. However, for T > there is generally 
sufficiently many steps between restarts for the history to grow to full length, 
and the rate of convergence does improve somewhat. 

For the Z4-Z3 model the energy difference between the highest occupied 
and lowest unoccupied orbital is 1.69 x 10~ , see Figure [3] and Table [a] This 
difference is comparatively small and could explain the poor performance of the 
QN method, particularly for T = 0. In Figure |8] the convergence rate of the 
optimization procedures for Z4-Z3 for T — 0.3, 0.5, 0.7, and the convergence rate 
for T = is included for reference. At T = 0.3 the rate of convergence for the 
QN method is considerably improved and the convergence rate remains superior 
to r = for T = 0.5 and T = 0.7. The elevated temperature broadens the Fermi 
surface and this could explain the improved convergence at T = 0.3, while the 
convergence of higher energy orbitals makes the problem more challenging at 
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Figure 4: Energy convergence for at varying temperatures. 
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Figure 5: Energy convergence for ^3-^2 at varying temperatures. 
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Figure 6: Energy convergence for Za^-Zj, at varying temperatures. 
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Figure 7: Energy convergence of the simultaneous QN method with restarts for 
Z4-Z3 at varying temperatures. 
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Figure 8: Energy convergence for Z^-Z-^ at varying temperatures T < 1. 
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higher temperatures. This would also explain the decreasing performance of 
NLCG for higher temperatures. 



Table 1: Orbital energy levels with occupation for Z2 at varying temperatures. 
The same data is graphically presented in Figure [T| 
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Table 2: Orbital energy levels with occupation for Z^-Z2 at varying tempera- 
tures with orbital energies shifted by +5. The same data is graphically presented 

in F igure [2| 
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4 Conclusion 

We have presented two schemes for energy optimization of ensemble DFT com- 
putations. The updates take the problem constraints into account and permits 
us to use information obtained from previous evaluations of the target functional 
and gradients to improve rate of convergence. We have further demonstrated 
the methods numerically on a model problem inspired by the electronic struc- 
ture theory and compared simultaneous and sequential schemes based on the 
QN and NLCG methods. 
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The ensemble model successfully concentrates occupation to low energy or- 
bitals at low temperatures, and gradually increases occupation of higher energy 
orbitals at increasing temperature to increase the entropy of the system. Opti- 
mization of the occupation numbers also enables ensemble DFT calculations to 
automatically handle degenerate and near degenerate orbitals at T = 0, which 
are challenging for methods that construct the electron density by the Auf- 
bau principle. Furthermore, is seems possible to broaden the Fermi surface by 
increasing temperature to accelerate convergence of small gap systems. 

Simultaneous optimization schemes provide improved convergence compared 
to sequential approaches for both the NLCG and QN methods. While the NLCG 
and QN methods are often comparable in performance, the NLCG method is 
overall more robust. In contrast. Reference found that QN method is more 
robust than the NLCG method. It is possible that the quadratic approximate 
line search gives a better result for the model problem. As the NLCG method 
depends heavily on a high quality line search this might provide a possible 
explanation. In the present case, the QN method performs poorly for problems 
with frequent restarts and while this effect does not fully explain the lack of 
convergence it can be used as a problem indicator. 



Table 3: Orbital energy levels with occupation for Z4-Z3 at varying tempera- 
tures with orbital energies shifted by -1-10. The same data is graphically pre- 

sent ed in Figure [3} 

E Occ. (T=0) Occ. (T=l) Occ. (T=2) Occ. (T=3) 



4.327524 


1, 


.000000 


1, 


.000000 


1, 


.000000 


1, 


.000000 


17.873544 


1, 


.000000 


1, 


.000000 


1, 


.000000 


1, 


.000000 


25.639021 


1, 


.000000 


1, 


.000000 


1, 


.000000 


1, 


.000000 


38.541992 


1, 


.000000 


1, 


.000000 


1. 


.000000 


1, 


.000000 


47.960214 


1, 


.000000 


1, 


.000000 


0, 


.999983 


0, 


.994464 


48.278074 


1, 


.000000 


1, 


.000000 


0. 


.999841 


0, 


.993917 


63.300484 


1 


.000000 





.669980 


0. 


.591678 


0, 


.564432 


64.759294 


0, 


.000000 


0, 


.330020 


0, 


.408498 


0, 


.439697 


78.938961 


0, 


.000000 


0, 


.000000 


0, 


.000000 


0, 


.005937 


82.626842 


0, 


.000000 


0, 


.000000 


0, 


.000000 


0, 


.001553 


89.029063 


0, 


.000000 


0, 


.000000 


0, 


.000000 


0, 


.000000 


95.252731 


0, 


.000000 


0, 


.000000 


0, 


.000000 


0, 


.000000 


98.574017 


0, 


.000000 


0, 


.000000 


0, 


.000000 


0, 


.000000 
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